In this report, we reproduce the analyses in the follow-up behavioral study 2.

prep data

First, we load the relevant packages, define functions and plotting aesthetics, and load and tidy the data.

load packages

if(!require('pacman')) {
    install.packages('pacman')
}

pacman::p_load(tidyverse, purrr, fs, knitr, lmerTest, ggeffects, kableExtra, boot, devtools, EMAtools, install = TRUE)
devtools::install_github("hadley/emo")

define functions

# MLM results table function
table_model = function(model_data, eff_size = TRUE, word_count = TRUE) {
  
  results = model_data %>%
    broom.mixed::tidy(conf.int = TRUE) %>%
    filter(effect == "fixed") %>%
    rename("SE" = std.error,
           "t" = statistic,
           "p" = p.value) %>%
    select(-group, -effect) %>%
    mutate_at(vars(-contains("term"), -contains("p")), round, 2) %>%
    mutate(term = gsub("article_cond", "", term),
           term = gsub("\\(Intercept\\)", "control", term),
           term = gsub("sharing_type", "sharing type", term),
           term = gsub("msg_rel_self_between", "self-relevance", term),
           term = gsub("msg_rel_social_between", "social relevance", term),
           term = gsub("grouptimed", "group (timed)", term),
           term = gsub("groupuntimed", "group (untimed)", term),
           term = gsub("contentclimate", "content (climate)", term),
           term = gsub("siteUSA", "sample (USA)", term),
           term = gsub("n_c", "word count", term),
           term = gsub(":", " x ", term),
           p = ifelse(p < .001, "< .001",
                      ifelse(p == 1, "1.000", gsub("0.(.*)", ".\\1", sprintf("%.3f", p)))),
           `b [95% CI]` = sprintf("%.2f [%0.2f, %.2f]", estimate, conf.low, conf.high)) 
  
  if (word_count == TRUE) {
    results = results %>%
      mutate(term = gsub("control", "intercept", term))
  }
  
  if (eff_size == TRUE) {
    eff_size = lme.dscore(model_data, data = data, type = "lme4") %>%
      rownames_to_column(var = "term") %>%
      mutate(term = gsub("article_cond", "", term),
             term = gsub("article_cond", "", term),
             term = gsub("\\(Intercept\\)", "control", term),
             term = gsub("sharing_type", "sharing type", term),
             term = gsub("msg_rel_self_between", "self-relevance", term),
             term = gsub("msg_rel_social_between", "social relevance", term),
             term = gsub("contentclimate", "content (climate)", term),
             term = gsub(":", " x ", term),
             d = sprintf("%.2f", d)) %>%
      select(term, d)
    
    results %>%
      left_join(., eff_size) %>%
      mutate(d = ifelse(is.na(d), "--", d)) %>%
      select(term, `b [95% CI]`, d, df, t, p) %>%
      kable() %>%
      kableExtra::kable_styling()
    
  } else {
    results %>%
      select(term, `b [95% CI]`, df, t, p) %>%
      kable() %>%
      kableExtra::kable_styling()
  }
}

# simple effects function
simple_effects = function(model, sharing = FALSE) {
  if(sharing == FALSE) {
    results = emmeans::contrast(emmeans::emmeans(model, ~ article_cond | group),
                            "revpairwise", by = "group", adjust = "none") %>%
      data.frame() %>%
      filter(grepl("control", contrast)) %>%
      select(contrast, group, estimate, p.value)
  } else {
    results = emmeans::contrast(emmeans::emmeans(model, ~ article_cond | group + sharing_type),
                            "revpairwise", by = "group", adjust = "none") %>%
      data.frame() %>%
      filter(grepl("- control", contrast)) %>%
      filter(!grepl("^control", contrast)) %>%
      extract(contrast, c("exp_sharing", "control_sharing"), ".* (0|1) - control (0|1)", remove = FALSE) %>%
      filter(exp_sharing == control_sharing) %>%
      mutate(sharing_type = ifelse(exp_sharing == 0, "broadcast", "narrowcast"),
             contrast = gsub("0|1", "", contrast)) %>%
      select(contrast, sharing_type, group, estimate, p.value)
  }
  
  results %>%
    mutate(p.value = ifelse(p.value < .001, "< .001",
                      ifelse(p.value == 1, "1.000", gsub("0.(.*)", ".\\1", sprintf("%.3f", p.value))))) %>%
    kable(digits = 2) %>%
    kableExtra::kable_styling()
}

define aesthetics

palette_condition = c("#ee9b00", "#bb3e03", "#005f73")
palette_dv = c("#ee9b00", "#005f73", "#56282D")
palette_sharing = c("#0a9396", "#ee9b00")

plot_aes = theme_minimal() +
  theme(legend.position = "top",
        legend.text = element_text(size = 12),
        text = element_text(size = 16, family = "Futura Medium"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.text = element_text(color = "black"),
        axis.line = element_line(colour = "black"),
        axis.ticks.y = element_blank())

load data

data = read.csv("../data/study2_data.csv", stringsAsFactors = FALSE) %>%
  mutate(article_cond = ifelse(article_cond == "social", "other", article_cond))

n_words = read.csv("../data/study2_n_words.csv", stringsAsFactors = FALSE) %>%
  mutate(article_cond = ifelse(article_cond == "social", "other", article_cond))

descriptives

group ns

Sample size by group

data %>%
  select(group, SID) %>%
  unique() %>%
  group_by(group) %>%
  summarize(n = n()) %>%
  kable() %>%
  kable_styling()
group n
comment 131
timed 159
untimed 169

ratings

Summarize means and SDs

data %>%
  gather(variable, value, msg_share, msg_rel_self, msg_rel_social) %>%
  group_by(variable) %>%
  summarize(M = mean(value, na.rm = TRUE),
            SD = sd(value, na.rm = TRUE)) %>%
  mutate(variable = ifelse(variable == "msg_rel_self", "self-relevance",
                    ifelse(variable == "msg_rel_social", "social relevance", "sharing intention"))) %>%
  kable(digits = 2) %>%
  kableExtra::kable_styling()
variable M SD
self-relevance 49.40 34.17
social relevance 54.54 33.03
sharing intention 31.65 32.94

manipulation checks

H2

Do the manipulations increase relevance? Is this effect stronger in the comment group?

self-relevance

✅ H2a: Self-focused intervention (compared to control) will increase self-relevance

We replicate our previous work in the comment group: the self-focused condition increases self-relevance compared to the control

✅ This effect is smaller for both the timed and untimed groups

mod_h2a = lmer(msg_rel_self ~ 1 + article_cond * group + (1 + article_cond | SID),
              data = filter(data, sharing_type == 1),
              control = lmerControl(optimizer = "bobyqa"))

model table

table_model(mod_h2a, eff_size = FALSE)
term b [95% CI] df t p
intercept 49.62 [45.38, 53.86] 451.69 23.01 < .001
other 4.05 [0.54, 7.57] 447.10 2.26 .024
self 14.61 [10.70, 18.53] 445.68 7.34 < .001
group (timed) -3.13 [-8.81, 2.54] 449.32 -1.09 .278
group (untimed) -4.53 [-10.13, 1.07] 450.27 -1.59 .113
other x group (timed) -5.71 [-10.41, -1.00] 445.43 -2.38 .018
self x group (timed) -9.94 [-15.16, -4.71] 445.32 -3.74 < .001
other x group (untimed) -3.33 [-7.98, 1.31] 445.67 -1.41 .159
self x group (untimed) -11.88 [-17.04, -6.72] 445.70 -4.52 < .001

simple effects by group

simple_effects(mod_h2a)
contrast group estimate p.value
other - control comment 4.05 .024
self - control comment 14.61 < .001
other - control timed -1.65 .298
self - control timed 4.68 .008
other - control untimed 0.72 .642
self - control untimed 2.73 .111

summary

summary(mod_h2a)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: msg_rel_self ~ 1 + article_cond * group + (1 + article_cond |  
##     SID)
##    Data: filter(data, sharing_type == 1)
## Control: lmerControl(optimizer = "bobyqa")
## 
## REML criterion at convergence: 51781.6
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.91268 -0.68072  0.04976  0.68791  2.92129 
## 
## Random effects:
##  Groups   Name              Variance Std.Dev. Corr       
##  SID      (Intercept)       394.93   19.873              
##           article_condother  28.66    5.353   -0.12      
##           article_condself  121.38   11.017   -0.22  0.49
##  Residual                   735.56   27.121              
## Number of obs: 5378, groups:  SID, 454
## 
## Fixed effects:
##                                Estimate Std. Error      df t value
## (Intercept)                      49.619      2.156 451.693  23.014
## article_condother                 4.054      1.791 447.103   2.264
## article_condself                 14.614      1.991 445.684   7.340
## grouptimed                       -3.135      2.886 449.319  -1.086
## groupuntimed                     -4.528      2.849 450.270  -1.589
## article_condother:grouptimed     -5.707      2.394 445.433  -2.384
## article_condself:grouptimed      -9.936      2.660 445.322  -3.735
## article_condother:groupuntimed   -3.335      2.365 445.668  -1.410
## article_condself:groupuntimed   -11.879      2.627 445.702  -4.521
##                                            Pr(>|t|)    
## (Intercept)                    < 0.0000000000000002 ***
## article_condother                          0.024040 *  
## article_condself                   0.00000000000102 ***
## grouptimed                                 0.277973    
## groupuntimed                               0.112677    
## article_condother:grouptimed               0.017544 *  
## article_condself:grouptimed                0.000212 ***
## article_condother:groupuntimed             0.159199    
## article_condself:groupuntimed      0.00000789291645 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##                   (Intr) artcl_cndt artcl_cnds grptmd grpntm artcl_cndthr:grpt
## artcl_cndth       -0.421                                                      
## artcl_cndsl       -0.440  0.488                                               
## grouptimed        -0.747  0.314      0.329                                    
## groupuntimd       -0.757  0.318      0.333      0.565                         
## artcl_cndthr:grpt  0.315 -0.748     -0.365     -0.416 -0.238                  
## artcl_cndslf:grpt  0.329 -0.365     -0.748     -0.439 -0.249  0.486           
## artcl_cndthr:grpn  0.319 -0.757     -0.369     -0.238 -0.417  0.566           
## artcl_cndslf:grpn  0.333 -0.369     -0.758     -0.249 -0.440  0.276           
##                   artcl_cndslf:grpt artcl_cndthr:grpn
## artcl_cndth                                          
## artcl_cndsl                                          
## grouptimed                                           
## groupuntimd                                          
## artcl_cndthr:grpt                                    
## artcl_cndslf:grpt                                    
## artcl_cndthr:grpn  0.276                             
## artcl_cndslf:grpn  0.567             0.486

social relevance

✅ H2b: Other-focused intervention (compared to control) will increase social relevance

We replicate our previous work in the comment group: the other-focused condition increases social relevance compared to the control

✅ This effect is smaller for both the timed and untimed groups

mod_h2b = lmer(msg_rel_social ~ 1 + article_cond * group + (1 + article_cond | SID),
              data = filter(data, sharing_type == 1),
              control = lmerControl(optimizer = "bobyqa"))

model table

table_model(mod_h2b, eff_size = FALSE)
term b [95% CI] df t p
intercept 53.33 [49.01, 57.64] 453.40 24.28 < .001
other 14.90 [11.41, 18.38] 447.67 8.40 < .001
self 10.49 [6.93, 14.04] 444.92 5.80 < .001
group (timed) -2.62 [-8.40, 3.16] 450.98 -0.89 .373
group (untimed) -4.06 [-9.76, 1.65] 452.03 -1.40 .163
other x group (timed) -14.43 [-19.09, -9.77] 445.92 -6.09 < .001
self x group (timed) -7.44 [-12.19, -2.70] 444.54 -3.08 .002
other x group (untimed) -10.72 [-15.32, -6.12] 446.17 -4.58 < .001
self x group (untimed) -8.05 [-12.75, -3.36] 444.82 -3.38 < .001

simple effects by group

simple_effects(mod_h2b)
contrast group estimate p.value
other - control comment 14.90 < .001
self - control comment 10.49 < .001
other - control timed 0.47 .765
self - control timed 3.04 .057
other - control untimed 4.18 .006
self - control untimed 2.43 .118

summary

summary(mod_h2b)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: msg_rel_social ~ 1 + article_cond * group + (1 + article_cond |  
##     SID)
##    Data: filter(data, sharing_type == 1)
## Control: lmerControl(optimizer = "bobyqa")
## 
## REML criterion at convergence: 50967.6
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.12951 -0.57636  0.08024  0.63937  3.08780 
## 
## Random effects:
##  Groups   Name              Variance Std.Dev. Corr       
##  SID      (Intercept)       448.17   21.170              
##           article_condother  82.70    9.094   -0.16      
##           article_condself   97.11    9.855   -0.29  0.49
##  Residual                   612.31   24.745              
## Number of obs: 5378, groups:  SID, 454
## 
## Fixed effects:
##                                Estimate Std. Error      df t value
## (Intercept)                      53.326      2.196 453.400  24.285
## article_condother                14.897      1.773 447.672   8.403
## article_condself                 10.487      1.808 444.918   5.799
## grouptimed                       -2.621      2.940 450.982  -0.892
## groupuntimed                     -4.058      2.902 452.033  -1.398
## article_condother:grouptimed    -14.427      2.371 445.916  -6.086
## article_condself:grouptimed      -7.444      2.416 444.536  -3.081
## article_condother:groupuntimed  -10.718      2.342 446.171  -4.577
## article_condself:groupuntimed    -8.055      2.386 444.825  -3.375
##                                            Pr(>|t|)    
## (Intercept)                    < 0.0000000000000002 ***
## article_condother              0.000000000000000584 ***
## article_condself               0.000000012657137565 ***
## grouptimed                                 0.373035    
## groupuntimed                               0.162702    
## article_condother:grouptimed   0.000000002499743111 ***
## article_condself:grouptimed                0.002189 ** 
## article_condother:groupuntimed 0.000006126021397687 ***
## article_condself:groupuntimed              0.000802 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##                   (Intr) artcl_cndt artcl_cnds grptmd grpntm artcl_cndthr:grpt
## artcl_cndth       -0.393                                                      
## artcl_cndsl       -0.436  0.501                                               
## grouptimed        -0.747  0.293      0.326                                    
## groupuntimd       -0.757  0.297      0.330      0.565                         
## artcl_cndthr:grpt  0.294 -0.748     -0.375     -0.388 -0.222                  
## artcl_cndslf:grpt  0.327 -0.375     -0.749     -0.436 -0.247  0.500           
## artcl_cndthr:grpn  0.297 -0.757     -0.380     -0.222 -0.389  0.566           
## artcl_cndslf:grpn  0.331 -0.380     -0.758     -0.247 -0.437  0.284           
##                   artcl_cndslf:grpt artcl_cndthr:grpn
## artcl_cndth                                          
## artcl_cndsl                                          
## grouptimed                                           
## groupuntimd                                          
## artcl_cndthr:grpt                                    
## artcl_cndslf:grpt                                    
## artcl_cndthr:grpn  0.284                             
## artcl_cndslf:grpn  0.567             0.500

combined plot

# generate predicted values
predicted_h2a = ggeffects::ggpredict(mod_h2a, c("article_cond", "group")) %>%
              data.frame() %>%
  mutate(model = "self-relevance")

predicted_h2b = ggeffects::ggpredict(mod_h2b, c("article_cond", "group")) %>%
              data.frame() %>%
  mutate(model = "social relevance")

# manipulation check plot
bind_rows(predicted_h2a, predicted_h2b) %>%
  mutate(x = factor(x, levels = c("self", "control", "other"))) %>%
  ggplot(aes(x = group, y = predicted, color = x)) +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high), position = position_dodge(.5), size = 1) +
  facet_grid(~ model) +
  coord_flip() +
  scale_color_manual(name = "", values = palette_condition) +
  labs(x = "", y = "\nmean predicted relevance rating") +
  plot_aes +
  theme(legend.position = "top")

H5

Do the manipulations increase sharing intentions? Is this effect stronger in the comment group?

narrowcasting only

Here we focus on narrowcasting only since that is the type of sharing we used in fMRI study 1.

✅ H5a: Self-focused intervention (compared to control) will increase sharing intentions

✅ H5b: Other-focused intervention (compared to control) will increase sharing intentions

We replicate our previous work in the comment group: the self- and other-focused conditions increase sharing intentions compared to the control, and these effects are stronger for narrowcast compared to broadcasting sharing intentions

✅ These effects were smaller for both the timed and untimed groups

mod_h5 = lmer(msg_share ~ 1 + article_cond*group + (1 + article_cond | SID),
              data = filter(data, sharing_type == 1),
              control = lmerControl(optimizer = "bobyqa"))

plot

# generate predicted values
predicted_h5 = ggeffects::ggpredict(mod_h5, c("article_cond", "group")) %>%
              data.frame() %>%
  mutate(model = "sharing")

# causal analysis plot
predicted_h5 %>%
  mutate(x = factor(x, levels = c("self", "control", "other"))) %>%
  ggplot(aes(x = group, y = predicted, color = x)) +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high), position = position_dodge(.5), size = 1) +
  coord_flip() +
  scale_color_manual(name = "", values = palette_condition) +
  labs(x = "", y = "\nmean predicted sharing intention rating") +
  plot_aes +
  theme(legend.position = "top")

model table

table_model(mod_h5, eff_size = FALSE)
term b [95% CI] df t p
intercept 32.44 [27.76, 37.13] 453.83 13.62 < .001
other 14.74 [11.40, 18.07] 448.37 8.69 < .001
self 9.54 [6.31, 12.78] 445.40 5.80 < .001
group (timed) 1.10 [-5.18, 7.37] 451.38 0.34 .731
group (untimed) 0.01 [-6.18, 6.20] 452.53 0.00 .997
other x group (timed) -13.78 [-18.24, -9.32] 446.43 -6.07 < .001
self x group (timed) -8.49 [-12.81, -4.17] 445.02 -3.86 < .001
other x group (untimed) -13.46 [-17.87, -9.05] 446.74 -6.00 < .001
self x group (untimed) -9.37 [-13.64, -5.10] 445.28 -4.32 < .001

simple effects by group

simple_effects(mod_h5, sharing = FALSE)
contrast group estimate p.value
other - control comment 14.74 < .001
self - control comment 9.54 < .001
other - control timed 0.96 .524
self - control timed 1.05 .471
other - control untimed 1.28 .382
self - control untimed 0.17 .904

summary

summary(mod_h5)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: msg_share ~ 1 + article_cond * group + (1 + article_cond | SID)
##    Data: filter(data, sharing_type == 1)
## Control: lmerControl(optimizer = "bobyqa")
## 
## REML criterion at convergence: 50202.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.4366 -0.5035 -0.0986  0.5111  3.7779 
## 
## Random effects:
##  Groups   Name              Variance Std.Dev. Corr       
##  SID      (Intercept)       585.72   24.202              
##           article_condother 106.11   10.301   -0.01      
##           article_condself   83.45    9.135   -0.09  0.31
##  Residual                   501.03   22.384              
## Number of obs: 5378, groups:  SID, 454
## 
## Fixed effects:
##                                 Estimate Std. Error        df t value
## (Intercept)                     32.44303    2.38270 453.82785  13.616
## article_condother               14.73944    1.69708 448.37035   8.685
## article_condself                 9.54382    1.64571 445.40060   5.799
## grouptimed                       1.09804    3.19211 451.38037   0.344
## groupuntimed                     0.01214    3.14968 452.53259   0.004
## article_condother:grouptimed   -13.78025    2.26970 446.43299  -6.071
## article_condself:grouptimed     -8.49349    2.19857 445.01508  -3.863
## article_condother:groupuntimed -13.45896    2.24205 446.74008  -6.003
## article_condself:groupuntimed   -9.37331    2.17184 445.28231  -4.316
##                                            Pr(>|t|)    
## (Intercept)                    < 0.0000000000000002 ***
## article_condother              < 0.0000000000000002 ***
## article_condself                      0.00000001264 ***
## grouptimed                                 0.731018    
## groupuntimed                               0.996925    
## article_condother:grouptimed          0.00000000271 ***
## article_condself:grouptimed                0.000129 ***
## article_condother:groupuntimed        0.00000000401 ***
## article_condself:groupuntimed         0.00001960750 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##                   (Intr) artcl_cndt artcl_cnds grptmd grpntm artcl_cndthr:grpt
## artcl_cndth       -0.265                                                      
## artcl_cndsl       -0.301  0.451                                               
## grouptimed        -0.746  0.198      0.225                                    
## groupuntimd       -0.756  0.201      0.228      0.565                         
## artcl_cndthr:grpt  0.198 -0.748     -0.337     -0.260 -0.150                  
## artcl_cndslf:grpt  0.226 -0.338     -0.749     -0.300 -0.171  0.450           
## artcl_cndthr:grpn  0.201 -0.757     -0.341     -0.150 -0.261  0.566           
## artcl_cndslf:grpn  0.228 -0.342     -0.758     -0.170 -0.302  0.256           
##                   artcl_cndslf:grpt artcl_cndthr:grpn
## artcl_cndth                                          
## artcl_cndsl                                          
## grouptimed                                           
## groupuntimd                                          
## artcl_cndthr:grpt                                    
## artcl_cndslf:grpt                                    
## artcl_cndthr:grpn  0.256                             
## artcl_cndslf:grpn  0.567             0.450

condition effects by sharing type

Here we include both narrowcast and broadcast sharing, and assess potential interactions.

✅ H5a: Self-focused intervention (compared to control) will increase sharing intentions

✅ H5b: Other-focused intervention (compared to control) will increase sharing intentions

We replicate our previous work in the comment group: the self and social conditions increase sharing intentions compared to the control, and these effects are stronger for narrowcast compared to broadcasting sharing intentions

✅ These effects were smaller for both the timed and untimed groups

mod_h5_sharing = lmer(msg_share ~ 1 + article_cond*sharing_type*group + (1 + sharing_type | SID),
              data = data,
              control = lmerControl(optimizer = "bobyqa"))

plot

# generate predicted values
predicted_h5_sharing = ggeffects::ggpredict(mod_h5_sharing, c("article_cond", "sharing_type", "group")) %>%
              data.frame() %>%
  mutate(group = ifelse(group == "0", "broadcast sharing", "narrowcast sharing"),
         facet = ifelse(grepl("time", facet), sprintf("reflect:\n%s", facet), "comment"),
         facet = factor(facet, levels = c("reflect:\ntimed", "reflect:\nuntimed", "comment")))

# causal analysis plot
predicted_h5_sharing %>%
  mutate(group = gsub(" sharing", "", group)) %>%
  ggplot(aes(x = facet, y = predicted, color = x)) +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high), position = position_dodge(.5), size = 1) +
  facet_grid(~ group) +
  coord_flip() +
  scale_color_manual(name = "", values = palette_condition) +
  labs(x = "", y = "\nmean predicted sharing intention rating") +
  plot_aes +
  theme(legend.position = "top")

model table

table_model(mod_h5_sharing, eff_size = FALSE)
term b [95% CI] df t p
intercept 24.65 [19.85, 29.45] 559.85 10.09 < .001
other 6.47 [3.87, 9.08] 9875.19 4.87 < .001
self 5.42 [2.81, 8.02] 9846.67 4.07 < .001
sharing type 7.78 [4.00, 11.57] 973.04 4.03 < .001
group (timed) 1.87 [-4.56, 8.30] 556.66 0.57 .568
group (untimed) 1.60 [-4.74, 7.95] 558.28 0.50 .620
other x sharing type 8.29 [4.61, 11.97] 9892.38 4.42 < .001
self x sharing type 4.14 [0.45, 7.82] 9851.42 2.20 .028
other x group (timed) -4.80 [-8.29, -1.32] 9860.03 -2.71 .007
self x group (timed) -5.25 [-8.73, -1.77] 9848.22 -2.95 .003
other x group (untimed) -5.33 [-8.77, -1.89] 9861.83 -3.04 .002
self x group (untimed) -4.29 [-7.73, -0.85] 9852.85 -2.45 .014
sharing type x group (timed) -0.75 [-5.81, 4.32] 963.48 -0.29 .773
sharing type x group (untimed) -1.63 [-6.64, 3.37] 967.36 -0.64 .521
other x sharing type x group (timed) -9.00 [-13.93, -4.08] 9869.82 -3.59 < .001
self x sharing type x group (timed) -3.26 [-8.19, 1.66] 9853.18 -1.30 .194
other x sharing type x group (untimed) -8.12 [-12.99, -3.26] 9871.90 -3.28 .001
self x sharing type x group (untimed) -5.07 [-9.93, -0.20] 9860.72 -2.04 .041

simple effects by group

simple_effects(mod_h5_sharing, sharing = TRUE)
contrast sharing_type group estimate p.value
other - control broadcast comment 6.47 < .001
self - control broadcast comment 5.42 < .001
other - control narrowcast comment 14.77 < .001
self - control narrowcast comment 9.55 < .001
other - control broadcast timed 1.67 .156
self - control broadcast timed 0.17 .887
other - control narrowcast timed 0.96 .416
self - control narrowcast timed 1.04 .378
other - control broadcast untimed 1.14 .318
self - control broadcast untimed 1.12 .326
other - control narrowcast untimed 1.31 .252
self - control narrowcast untimed 0.19 .868

summary

summary(mod_h5_sharing)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## msg_share ~ 1 + article_cond * sharing_type * group + (1 + sharing_type |  
##     SID)
##    Data: data
## Control: lmerControl(optimizer = "bobyqa")
## 
## REML criterion at convergence: 98063.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.0263 -0.4780 -0.0671  0.4288  4.3371 
## 
## Random effects:
##  Groups   Name         Variance Std.Dev. Corr 
##  SID      (Intercept)  641.6    25.33         
##           sharing_type 243.5    15.60    -0.35
##  Residual              436.3    20.89         
## Number of obs: 10756, groups:  SID, 454
## 
## Fixed effects:
##                                              Estimate Std. Error        df
## (Intercept)                                   24.6488     2.4422  559.8511
## article_condother                              6.4747     1.3285 9875.1864
## article_condself                               5.4175     1.3302 9846.6746
## sharing_type                                   7.7824     1.9297  973.0365
## grouptimed                                     1.8716     3.2747  556.6600
## groupuntimed                                   1.6025     3.2297  558.2793
## article_condother:sharing_type                 8.2926     1.8770 9892.3783
## article_condself:sharing_type                  4.1360     1.8809 9851.4232
## article_condother:grouptimed                  -4.8044     1.7760 9860.0304
## article_condself:grouptimed                   -5.2506     1.7770 9848.2188
## article_condother:groupuntimed                -5.3299     1.7546 9861.8350
## article_condself:groupuntimed                 -4.2931     1.7554 9852.8500
## sharing_type:grouptimed                       -0.7462     2.5807  963.4809
## sharing_type:groupuntimed                     -1.6344     2.5484  967.3642
## article_condother:sharing_type:grouptimed     -9.0044     2.5103 9869.8216
## article_condself:sharing_type:grouptimed      -3.2645     2.5125 9853.1763
## article_condother:sharing_type:groupuntimed   -8.1247     2.4798 9871.8961
## article_condself:sharing_type:groupuntimed    -5.0696     2.4817 9860.7176
##                                             t value             Pr(>|t|)    
## (Intercept)                                  10.093 < 0.0000000000000002 ***
## article_condother                             4.874           0.00000111 ***
## article_condself                              4.073           0.00004681 ***
## sharing_type                                  4.033           0.00005942 ***
## grouptimed                                    0.572             0.567861    
## groupuntimed                                  0.496             0.619966    
## article_condother:sharing_type                4.418           0.00001006 ***
## article_condself:sharing_type                 2.199             0.027903 *  
## article_condother:grouptimed                 -2.705             0.006840 ** 
## article_condself:grouptimed                  -2.955             0.003136 ** 
## article_condother:groupuntimed               -3.038             0.002390 ** 
## article_condself:groupuntimed                -2.446             0.014475 *  
## sharing_type:grouptimed                      -0.289             0.772526    
## sharing_type:groupuntimed                    -0.641             0.521458    
## article_condother:sharing_type:grouptimed    -3.587             0.000336 ***
## article_condself:sharing_type:grouptimed     -1.299             0.193872    
## article_condother:sharing_type:groupuntimed  -3.276             0.001055 ** 
## article_condself:sharing_type:groupuntimed   -2.043             0.041103 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

combined plot

bind_rows(predicted_h2a, predicted_h2b, predicted_h5) %>%
  mutate(model = factor(model, levels = c("self-relevance", "social relevance", "sharing")),
         x = factor(x, levels = c("self", "control", "other")),
         group = ifelse(group == "timed", "reflect:\ntimed",
                 ifelse(group == "untimed", "reflect:\nuntimed", "comment")),
         group = factor(group, levels = c("reflect:\ntimed", "reflect:\nuntimed", "comment"))) %>%
  ggplot(aes(x = group, y = predicted, color = x)) +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high), position = position_dodge(.5), size = 1.5) +
  facet_grid(~ model) +
  coord_flip() +
  scale_color_manual(name = "", values = palette_condition) +
  labs(x = "", y = "\npredicted rating") +
  plot_aes +
  theme(legend.position = "top")

word count analyses

Test whether word count is higher in the experimental conditions, and whether it is positively associated with self and social relevance, and sharing intention ratings.

descriptives

words_ratings = n_words %>%
  left_join(., data) %>%
  ungroup() %>%
  mutate(n_c = n - mean(n, na.rm = TRUE))

n_words %>%
  group_by(article_cond) %>%
  summarize(mean = mean(n, na.rm = TRUE),
            sd = sd(n, na.rm = TRUE),
            min = min(n, na.rm = TRUE),
            max = max(n, na.rm = TRUE)) %>%
  kable(digits = 2) %>%
  kableExtra::kable_styling()
article_cond mean sd min max
control 13.77 7.33 3 72
other 17.17 9.32 3 69
self 18.14 10.43 3 72

condition effects

Is word count higher in the experimental conditions compared to the control condition?

✅ The word count is higher in the experimental conditions compared to the control condition

mod_words = lmer(n ~ 1 + article_cond + (1 + article_cond | SID),
              data = n_words,
              control = lmerControl(optimizer = "bobyqa"))

plot

predicted = ggeffects::ggpredict(mod_words, c("article_cond")) %>%
              data.frame()

predicted %>%
  ggplot(aes(x = "", y = predicted, color = x)) +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high), position = position_dodge(.5), size = 1) +
  coord_flip() +
  scale_color_manual(name = "", values = palette_condition) +
  labs(x = "", y = "\nmean predicted word count") +
  plot_aes +
  theme(legend.position = "top")

model table

table_model(mod_words, eff_size = FALSE)
term b [95% CI] df t p
intercept 13.76 [12.74, 14.79] 126.25 26.56 < .001
other 3.38 [2.33, 4.43] 124.53 6.38 < .001
self 4.33 [3.15, 5.52] 124.73 7.23 < .001

summary

summary(mod_words)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: n ~ 1 + article_cond + (1 + article_cond | SID)
##    Data: n_words
## Control: lmerControl(optimizer = "bobyqa")
## 
## REML criterion at convergence: 9994.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.9259 -0.5021 -0.1098  0.3879  6.4856 
## 
## Random effects:
##  Groups   Name              Variance Std.Dev. Corr     
##  SID      (Intercept)       24.92    4.992             
##           article_condother 17.94    4.236    0.20     
##           article_condself  27.51    5.245    0.33 0.71
##  Residual                   34.33    5.859             
## Number of obs: 1492, groups:  SID, 127
## 
## Fixed effects:
##                   Estimate Std. Error       df t value             Pr(>|t|)    
## (Intercept)        13.7628     0.5181 126.2507  26.563 < 0.0000000000000002 ***
## article_condother   3.3820     0.5302 124.5339   6.379      0.0000000031775 ***
## article_condself    4.3344     0.5993 124.7301   7.232      0.0000000000422 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) artcl_cndt
## artcl_cndth -0.142           
## artcl_cndsl -0.007  0.615

relevance

Is word count positively associated with self and social relevance ratings?

self-relevance

✅ Word count is positively associated with self-relevance ratings

mod_words_h1 = lmer(msg_rel_self ~ 1 + n_c + (1 + n_c | SID),
              data = filter(words_ratings, sharing_type == 1),
              control = lmerControl(optimizer = "bobyqa"))

plot

values = seq(-15, 60, 10)
predicted_self = ggeffects::ggpredict(mod_words_h1, terms = "n_c [values]") %>%
  data.frame()

predicted_self %>%
  ggplot(aes(x, predicted)) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .25) +
  geom_line(size = 1) +
  coord_cartesian(ylim = c(40, 110)) +
  labs(x = "\nword count (grand mean-centered)", y = "predicted self-relevance rating\n") +
  plot_aes

model table

table_model(mod_words_h1, eff_size = FALSE, word_count = TRUE)
term b [95% CI] df t p
intercept 55.77 [51.77, 59.78] 126.45 27.55 < .001
word count 0.49 [0.24, 0.74] 56.60 3.97 < .001

summary

summary(mod_words_h1)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: msg_rel_self ~ 1 + n_c + (1 + n_c | SID)
##    Data: filter(words_ratings, sharing_type == 1)
## Control: lmerControl(optimizer = "bobyqa")
## 
## REML criterion at convergence: 14325.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9287 -0.6177  0.1320  0.6665  2.7664 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr
##  SID      (Intercept) 446.8830 21.1396      
##           n_c           0.2678  0.5175  0.18
##  Residual             720.1887 26.8363      
## Number of obs: 1491, groups:  SID, 127
## 
## Fixed effects:
##             Estimate Std. Error       df t value             Pr(>|t|)    
## (Intercept)  55.7736     2.0248 126.4522  27.546 < 0.0000000000000002 ***
## n_c           0.4894     0.1232  56.6018   3.973             0.000203 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##     (Intr)
## n_c 0.106

social relevance

✅ Word count is positively associated with social relevance ratings

mod_words_h2 = lmer(msg_rel_social ~ 1 + n_c + (1 + n_c | SID),
              data = filter(words_ratings, sharing_type == 1),
              control = lmerControl(optimizer = "bobyqa"))

plot

values = seq(-15, 60, 10)
predicted_social = ggeffects::ggpredict(mod_words_h2, terms = "n_c [values]") %>%
  data.frame()

predicted_social %>%
  ggplot(aes(x, predicted)) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .25) +
  geom_line(size = 1) +
  coord_cartesian(ylim = c(40, 105)) +
  labs(x = "\nword count (grand mean-centered)", y = "predicted social relevance rating\n") +
  plot_aes

model table

table_model(mod_words_h2, eff_size = FALSE, word_count = TRUE)
term b [95% CI] df t p
intercept 61.98 [57.89, 66.07] 126.76 30.01 < .001
word count 0.46 [0.23, 0.69] 77.46 4.04 < .001

summary

summary(mod_words_h2)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: msg_rel_social ~ 1 + n_c + (1 + n_c | SID)
##    Data: filter(words_ratings, sharing_type == 1)
## Control: lmerControl(optimizer = "bobyqa")
## 
## REML criterion at convergence: 14008
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.0851 -0.4808  0.1363  0.6187  2.9247 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr
##  SID      (Intercept) 481.6772 21.9471      
##           n_c           0.2764  0.5257  0.21
##  Residual             566.7872 23.8073      
## Number of obs: 1491, groups:  SID, 127
## 
## Fixed effects:
##             Estimate Std. Error       df t value             Pr(>|t|)    
## (Intercept)  61.9800     2.0656 126.7621  30.005 < 0.0000000000000002 ***
## n_c           0.4629     0.1147  77.4581   4.038             0.000126 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##     (Intr)
## n_c 0.127

combined plot

data_raw = words_ratings %>%
  filter(sharing_type == 1) %>%
  select(SID, n_c, msg_rel_self, msg_rel_social) %>%
  gather(group, predicted, contains("msg")) %>%
  rename("x" = n_c) %>%
  mutate(group = ifelse(group == "msg_rel_self", "self","social"),
         group = factor(group, levels = c("self", "social")))
  
predicted_self %>%
  mutate(group = "self") %>%
  bind_rows(., predicted_social %>%  mutate(group = "social")) %>%
  mutate(group = factor(group, levels = c("self", "social"))) %>%
  ggplot(aes(x, predicted, color = group, fill = group)) +
  geom_point(data = data_raw, aes(x, predicted, color = group, fill = group), alpha = .25) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .25, color = NA) +
  geom_line(size = 2) +
  scale_x_continuous(breaks = seq(-10, 60, 10)) +
  scale_y_continuous(breaks = seq(0, 100, 25)) +
  scale_color_manual(values = c(palette_condition[1], palette_condition[3]), name = "") + 
  scale_fill_manual(values = c(palette_condition[1], palette_condition[3]), name = "") + 
  labs(x = "\nword count (grand mean-centered)", y = "predicted relevance rating\n") +
  plot_aes +
  theme(legend.position = c(.85, .21))

sharing

Is word count positively associated with sharing intention ratings?

narrowcasting only

Here we focus on narrowcasting only since that is the type of sharing we used in fMRI study 1.

✅ Word count is positively associated with narrowcast sharing intentions

mod_words_h3 = lmer(msg_share ~ 1 + n_c + (1 + n_c | SID),
              data = filter(words_ratings, sharing_type == 1),
              control = lmerControl(optimizer = "bobyqa"))

plot

values = seq(-15, 60, 10)
predicted_sharing = ggeffects::ggpredict(mod_words_h3, terms = c("n_c [values]")) %>%
  data.frame()

predicted_sharing %>%
  ggplot(aes(x, predicted)) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .25, color = NA) +
  geom_line(size = 1) +
  scale_color_manual(values = palette_sharing, name = "") +
  scale_fill_manual(values = palette_sharing, name = "") +
  labs(x = "\nword count (grand mean-centered)", y = "predicted sharing intention rating\n") +
  plot_aes

model table

table_model(mod_words_h3, eff_size = FALSE, word_count = TRUE)
term b [95% CI] df t p
intercept 40.89 [36.25, 45.53] 126.13 17.43 < .001
word count 0.33 [0.07, 0.59] 52.69 2.56 .013

summary

summary(mod_words_h3)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: msg_share ~ 1 + n_c + (1 + n_c | SID)
##    Data: filter(words_ratings, sharing_type == 1)
## Control: lmerControl(optimizer = "bobyqa")
## 
## REML criterion at convergence: 14090
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.0304 -0.5457 -0.0914  0.5401  3.6751 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr
##  SID      (Intercept) 630.9952 25.1196      
##           n_c           0.5199  0.7211  0.09
##  Residual             581.2520 24.1092      
## Number of obs: 1491, groups:  SID, 127
## 
## Fixed effects:
##             Estimate Std. Error       df t value            Pr(>|t|)    
## (Intercept)  40.8898     2.3463 126.1264  17.427 <0.0000000000000002 ***
## n_c           0.3296     0.1285  52.6878   2.565              0.0132 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##     (Intr)
## n_c 0.093

sharing models by sharing type

✅ Word count is positively associated with sharing intentions (averaging across sharing types), but doesn’t differ by sharing type

mod_words_h4 = lmer(msg_share ~ 1 + n_c * sharing_type + (1 + n_c | SID),
              data = words_ratings,
              control = lmerControl(optimizer = "bobyqa"))

plot

values = seq(-20, 60, 10)
predicted = ggeffects::ggpredict(mod_words_h4, terms = c("n_c [values]", "sharing_type")) %>%
  data.frame() %>%
  mutate(group = ifelse(group == "0", "broadcast sharing", "narrowcast sharing"))

predicted %>%
  ggplot(aes(x, predicted, color = group, fill = group)) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .25, color = NA) +
  geom_line(size = 1) +
  scale_color_manual(values = palette_sharing, name = "") +
  scale_fill_manual(values = palette_sharing, name = "") +
  labs(x = "\nword count (grand mean-centered)", y = "predicted sharing intention rating\n") +
  plot_aes

model table

table_model(mod_words_h4, eff_size = FALSE, word_count = TRUE)
term b [95% CI] df t p
intercept 28.77 [24.23, 33.32] 135.12 12.53 < .001
word count 0.21 [0.01, 0.41] 86.66 2.05 .043
sharing type 12.09 [10.47, 13.70] 2719.42 14.68 < .001
word count x sharing type 0.07 [-0.10, 0.25] 2719.42 0.83 .408

summary

summary(mod_words_h4)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: msg_share ~ 1 + n_c * sharing_type + (1 + n_c | SID)
##    Data: words_ratings
## Control: lmerControl(optimizer = "bobyqa")
## 
## REML criterion at convergence: 27502.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1282 -0.5549 -0.0911  0.4457  3.3863 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr
##  SID      (Intercept) 616.4090 24.8276      
##           n_c           0.3193  0.5651  0.07
##  Residual             505.3583 22.4802      
## Number of obs: 2982, groups:  SID, 127
## 
## Fixed effects:
##                    Estimate Std. Error         df t value            Pr(>|t|)
## (Intercept)        28.77413    2.29676  135.11723  12.528 <0.0000000000000002
## n_c                 0.20936    0.10211   86.65880   2.050              0.0434
## sharing_type       12.08539    0.82333 2719.42381  14.679 <0.0000000000000002
## n_c:sharing_type    0.07330    0.08853 2719.42381   0.828              0.4077
##                     
## (Intercept)      ***
## n_c              *  
## sharing_type     ***
## n_c:sharing_type    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) n_c    shrng_
## n_c          0.069              
## sharing_typ -0.179  0.000       
## n_c:shrng_t  0.000 -0.434 -0.001

combined plot

data_raw = words_ratings %>%
  filter(sharing_type == 1) %>%
  select(SID, n_c, msg_rel_self, msg_rel_social, msg_share) %>%
  gather(group, predicted, contains("msg")) %>%
  rename("x" = n_c) %>%
  mutate(group = ifelse(group == "msg_rel_self", "self-relevance",
                 ifelse(group == "msg_rel_social", "social relevance", "sharing")),
         group = factor(group, levels = c("self-relevance", "social relevance", "sharing")))
  
predicted_self %>%
  mutate(group = "self-relevance") %>%
  bind_rows(., predicted_social %>%  mutate(group = "social relevance")) %>%
  bind_rows(., predicted_sharing %>%  mutate(group = "sharing")) %>%
  mutate(group = factor(group, levels = c("self-relevance", "social relevance", "sharing"))) %>%
  ggplot(aes(x, predicted, color = group, fill = group)) +
  geom_point(data = data_raw, aes(x, predicted, color = group, fill = group), alpha = .25) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .25, color = NA) +
  geom_line(size = 2) +
  scale_x_continuous(breaks = seq(-10, 60, 10)) +
  scale_y_continuous(breaks = seq(0, 100, 25)) +
  scale_color_manual(values = palette_dv, name = "") + 
  scale_fill_manual(values = palette_dv, name = "") + 
  labs(x = "\nword count (grand mean-centered)", y = "predicted rating\n") +
  plot_aes +
  theme(legend.position = c(.85, .21))

cite packages

report::cite_packages()
##   - Angelo Canty and Brian Ripley (2021). boot: Bootstrap R (S-Plus) Functions. R package version 1.3-28.
##   - Douglas Bates, Martin Maechler and Mikael Jagan (2023). Matrix: Sparse and Dense Matrix Classes and Methods. R package version 1.5-4. https://CRAN.R-project.org/package=Matrix
##   - Douglas Bates, Martin Maechler, Ben Bolker, Steve Walker (2015). Fitting Linear Mixed-Effects Models Using lme4. Journal of Statistical Software, 67(1), 1-48. doi:10.18637/jss.v067.i01.
##   - Evan Kleiman (2021). EMAtools: Data Management Tools for Real-Time Monitoring/Ecological Momentary Assessment Data. R package version 0.1.4. https://CRAN.R-project.org/package=EMAtools
##   - H. Wickham. ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York, 2016.
##   - Hadley Wickham (2019). stringr: Simple, Consistent Wrappers for Common String Operations. R package version 1.4.0. https://CRAN.R-project.org/package=stringr
##   - Hadley Wickham (2021). forcats: Tools for Working with Categorical Variables (Factors). R package version 0.5.1. https://CRAN.R-project.org/package=forcats
##   - Hadley Wickham and Maximilian Girlich (2022). tidyr: Tidy Messy Data. R package version 1.2.0. https://CRAN.R-project.org/package=tidyr
##   - Hadley Wickham, Jennifer Bryan and Malcolm Barrett (2021). usethis: Automate Package and Project Setup. R package version 2.1.5. https://CRAN.R-project.org/package=usethis
##   - Hadley Wickham, Jim Hester and Jennifer Bryan (2022). readr: Read Rectangular Text Data. R package version 2.1.2. https://CRAN.R-project.org/package=readr
##   - Hadley Wickham, Jim Hester, Winston Chang and Jennifer Bryan (2021). devtools: Tools to Make Developing R Packages Easier. R package version 2.4.3. https://CRAN.R-project.org/package=devtools
##   - Hadley Wickham, Romain François, Lionel Henry and Kirill Müller (2022). dplyr: A Grammar of Data Manipulation. R package version 1.0.9. https://CRAN.R-project.org/package=dplyr
##   - Hao Zhu (2021). kableExtra: Construct Complex Table with 'kable' and Pipe Syntax. R package version 1.3.4. https://CRAN.R-project.org/package=kableExtra
##   - Jim Hester, Hadley Wickham and Gábor Csárdi (2021). fs: Cross-Platform File System Operations Based on 'libuv'. R package version 1.5.2. https://CRAN.R-project.org/package=fs
##   - Kirill Müller and Hadley Wickham (2022). tibble: Simple Data Frames. R package version 3.1.8. https://CRAN.R-project.org/package=tibble
##   - Kuznetsova A, Brockhoff PB, Christensen RHB (2017). "lmerTest Package:Tests in Linear Mixed Effects Models." _Journal of StatisticalSoftware_, *82*(13), 1-26. doi: 10.18637/jss.v082.i13 (URL:https://doi.org/10.18637/jss.v082.i13).
##   - Lionel Henry and Hadley Wickham (2020). purrr: Functional Programming Tools. R package version 0.3.4. https://CRAN.R-project.org/package=purrr
##   - Lüdecke D (2018). "ggeffects: Tidy Data Frames of Marginal Effects fromRegression Models." _Journal of Open Source Software_, *3*(26), 772.doi: 10.21105/joss.00772 (URL: https://doi.org/10.21105/joss.00772).
##   - R Core Team (2021). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/.
##   - Rinker, T. W. & Kurkiewicz, D. (2017). pacman: Package Management for R. version 0.5.0. Buffalo, New York. http://github.com/trinker/pacman
##   - Wickham et al., (2019). Welcome to the tidyverse. Journal of Open Source Software, 4(43), 1686, https://doi.org/10.21105/joss.01686
##   - Yihui Xie (2021). knitr: A General-Purpose Package for Dynamic Report Generation in R. R package version 1.37.